Vortex lattice melting in layered superconductors with periodic columnar pins 
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The melting transition of the vortex lattice in highly anisotropic, layered superconductors with 
commensurate, periodic columnar pins is studied in a geometry where magnetic field and columnar 
pins are normal to the layers. Thermodynamic properties and equilibrium density distributions are 
obtained from numerical minimizations of an appropriate free-energy functional. We find a line of 
first-order transitions that ends at a critical point as the pin concentration is increased. A simple 
Landau theory providing a semi-quantitative explanation of the numerical results is proposed. 
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Thermodynamic and transport properties of type-II 
superconductors in the mixed phase can be altered in a 
controlled manner by the introduction of artificial pin- 
ning centers. Columnar pinning in high-temperature 
superconductors (HTSCs), produced by damage tracks 
from heavy-ion bombardment, has been extensively stud- 
ied experimentally and theoretically because 
such pinning enhances the critical current. A random ar- 
ray of such pins leads to a glassy phase, known as Bose 
glass P]. The interplay between the lattice constant of 
the pin array and the intervortex separation should pro- 
duce p| commensurability effects. Such effects have been 
observed in imaging experiments Q and magnetization 
measurements |Q on thin-film superconductors with pe- 
riodic arrays of artificially produced pinning centers. 

In this Letter, we investigate the effects of a commen- 
surate, periodic array of columnar pins on the vortex- 
lattice melting transition || in highly anisotropic, lay- 
ered, HTSCs in a magnetic field. Both the field and 
the columnar pins are assumed to be normal to the lay- 
ers. We consider values of the field for which the pin 
areal density is smaller than that of the vortex lines. Al- 
though vortex lattice melting in HTSCs in the absence 
of pinning has been studied extensively || , not much is 
known about the effects of columnar pins on this transi- 
tion. Since columnar defects produce strong pinning, at 
temperatures near the melting temperature of the vor- 
tex lattice in the pure system each defect should pin a 
vortex. However, interstitial vortices, present when there 
are more vortices than pins, should Q undergo a melting 
transition at a temperature slightly higher than the melt- 
ing point of the pure vortex lattice || . Evidence for such 
melting of interstitial vortices has been found in exper- 
iments [^JJtJ on thin-film superconductors with periodic 
pinning, but the thermodynamics of this transition has 
not been studied. 

Since the defect-pinned vortices produce a periodic po- 
tential for the interstitial ones, the melting transition 
of the latter provides a physical realization of three- 



dimensional melting in the presence of an external pe- 
riodic potential. For small concentrations of pinning 
centers, we find a first-order melting transition from a 
crystalline solid to an inhomogeneous liquid. As the pin 
concentration is increased, the transition temperature in- 
creases and the latent heat and the jump in the crys- 
talline order parameter at the transition decrease. This 
line of first-order transitions terminates at a critical point 
beyond which the thermodynamic transition is replaced 
by a sharp crossover. This critical point is a rare, ex- 
perimentally realizable example of continuous melting in 
three dimensions. We show that a simple Landau the- 
ory provides a semi-quantitative understanding of most 
of our results. Such melting transitions are of interest in 
other systems such as atoms adsorbed on crystalline sub- 
strates jicj, colloidal particles in interfering laser fields 
|]i"lf and arrays of optical traps Jl^] . Our work is of rele- 
vance to these systems also. 

Our study is based on a model free energy functional 
Jl3],[l4[ for a system of "pancake" vortices Q in a highly 
anisotropic layered superconductor. The commensurate 
array of columnar pins is accounted for by an appropri- 
ate "external potential" term Q in the free energy func- 
tional. By numerically minimizing a discretized form of 
this functional, we have studied the effects of periodic 
pinning on the structure and thermodynamics of the liq- 
uid and crystalline states of this system. We consider a 
layered superconductor with vanishingly small Josephson 
coupling between layers (vortices on different layers are 
coupled via their electromagnetic interaction only). In 
this limit of effectively infinite anisotropy, which is ap- 
propriate |l4| for extremely anisotropic Bi- and Tl-based 
compounds, the energy of a system of pancake vortices 
residing on the superconducting layers may be written as 
a sum of anisotropic two-body interactions v(n, r) where 
n is the layer separation and r is the separation in the 
plane of the layers. The Fourier transform of v is JTJ] 

= 2^rA 2 [fcj + (4/d 2 )sin 2 (fc z d/2)] 
1 1 k 2 ± [l + \ 2 k 2 ± + 4{\ 2 /d 2 ) sin 2 (k z d/ 2)}' { ' 
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with r = /3d$2 / / g7r 2 A 2 and p = i/ kBT _ Here, k z (k±) 
is the component of k perpendicular (parallel) to the 
layer plane, d is the layer spacing, A(T) the penetra- 
tion depth in the layer plane, and $o the flux quan- 
tum. The intralayer potential v(n = 0, r) is repulsive and 
oc ln(r), whereas the interlayer potential v(n ^ 0,r), also 
oc ln(r), is attractive, weaker than the intralayer poten- 
tial by the factor d/X, and decreases exponentially with 
n as e~ nd / x . We use parameters appropriate to BSCCO 
i.e. A(T = 0) = 1500A and d = 15A, and assume a 
two-fluid T dependence of A with T c (0) = 85K. 

In density functional theory |l5|-|l5|] , the free energy of 
a state in a density configuration specified by p(i, r), the 
time averaged areal density of vortices at point r on the 
ith layer, is given in terms of equilibrium correlation func- 
tions of the layered liquid of pancake vortices. We use the 
Ramakrishnan-Yussouff free energy functional |l3| ] which 
is known |l4j to provide a quantitatively correct descrip- 
tion of the melting transition in our system in the ab- 
sence of pinning. Since the potential produced by a set 
of straight columnar pins perpendicular to the layers is 
the same on every layer, the time-averaged density of 
vortices at any point r must be the same on all layers: 
p(i,r) = p(r) for all i. Then, the free energy per layer 
may be written in a two-dimensional form: 
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Here, 5p(r) = p(r) — po, F n is the free energy of the uni- 
form liquid of areal density po (= -B/$o where B is the 
magnetic induction), V p (r) is the pinning potential, and 
C( r ) = Sn^( n ' r )' wnere C(n,r) is the direct pair cor- 
relation function jl5| of a layered liquid of pancake vor- 
tices. We use the results for C(n, r) obtained jl4j from a 
hypernetted chain calculation ]l5[ ]. 

The pinning potential at point r is given by (3V p (r) = 
Vo{\r — Rj|), where the sum is over all pinning cen- 
ters located at the points {Rj} on a plane, and Vb(r), 
the potential at r due to a pinning center at the origin is 
assumed to have the truncated parabolic form 



V (r) = -aT(l-r 2 /r 2 ) 



(3) 



for r < ro and Vo(r) = if r > r$. Here, ro is the range 
and a is a "strength" parameter chosen to ensure that 
on the average a pinning center traps < 1 vortex in the 
temperature range of interest. 

We find the minima of the free energy of Eq. (Q) using 
a methodology quite similar to that in our earlier studies 
p6| , p7t of hard-sphere systems. We discretize space by 
defining density variables {pk} at the sites of a periodic 



grid, and use a gradient descent method |T(| to locate the 
minima of the free energy of Eq.(|^) written as a function 
of {pk}- The use of fast Fourier transforms in the calcu- 
lation of lattice sums pl| speeds the computations and 
allows studies of larger systems. To accommodate a tri- 
angular lattice, we use a triangular grid with periodic 
boundary conditions. 

We have performed our studies for B — 2kG and 3kG. 
We first checked the results for crystallization without 
pinning. For this purpose, the computational system was 
one triangular lattice unit cell with lattice constant a (all 
lengths are in units of ag, with ttoq = 1/po — $o/_B), 
and the spacing h of the computational grid was chosen 
to have the values a/n with n = 16, 32, 64 and 128. 
The free energies of the crystal obtained |Q for all these 
values of n are essentially the same, indicating that the 
effects of discretization are minimal for h < a/16. The 
equilibrium value a m of the lattice parameter a was de- 
termined by finding the value of a that minimizes the free 
energy at a given B and T. The value of a m is found to 
be slightly higher than the spacing of a perfect triangular 
lattice of density po. This reflects the well-known result 
H that the density of a vortex lattice increases slightly 
at melting. 

The free energy of the crystal crosses zero at the transi- 
tion temperature T c = 18.45K at B = 2kG. As expected, 
T c is slightly higher than that obtained from approximate 
treatments [ fTI] of the same free energy. The entropy 
change As per vortex is 0.29fcs, and the jump in the 
order parameter m, (m is the magnitude of the Fourier 
component of p(r) at the shortest reciprocal lattice vector 
of the triangular lattice) is Am = 0.52. The Lindemann 
parameter calculated from the density distribution |l9| ] 
at the crystalline minimum is C = 0.26 at melting. Very 
similar results were obtained for B — 3kG (T c = 15.10K, 
As = 0.28fc B , Am = 0.52, C = 0.25). The close agree- 
ment of these results with those of other studies [p|Jl4|l 
establishes the validity of our numerical approach. 

Next, we considered the effects of a single pin and a 
pair of pins on the liquid-state properties. The single-pin 
calculation was done to fix the value of the parameter a 
of Eq.(|). Using r = 0.1a (~ 60A for B =2 kG), we find 
that in the temperature range of interest (T w 15-25K), 
the integrated density inside the range of the pinning po- 
tential is close to one if a w 0.05-0.06 (V (r = 0) w 7-9). 
The density distribution near the center of the defect is 
gaussian, p(r) oc exp(— otYr 2 /r 2 ,), expected from Eq.([|). 
Since pinning of multiple vortices at radiation-induced 
defects is not found in experiments, higher values of a 
were not considered. Our two-pin calculations show |l9| ] 
the expected oscillatory behavior || of the free energy as 
a function of the distance between the pins. 

To study melting in the presence of a commensurate, 
periodic array of pins, we considered a triangular lattice 
of pins with spacing equal to la m where I is an integer. 
Thus the pin concentration is c = l/l 2 . The computa- 
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tional cell used was one unit cell of the pin lattice (which 
contains l 2 unit cells of the vortex lattice) with one pin 
located at the center of the first unit cell of the vor- 
tex lattice. The value of h was fixed at a m /64. The 
crystalline (or liquid) minimum of the free energy was 
located by starting the minimization from a crystalline 
(or liquid-like) initial state, usually that obtained at a 
nearby T. For both values of B and smaller c, (I = 10, 
8, 7 and 6) we found a first-order transition (two distinct 
local minima whose free energies cross at the transition 
temperature). Results for B — 2kG are shown in Fig|l]. 
As expected, the presence of columnar pins increases T c . 
The discontinuities in s and m decrease as c increases 
(inset in Fig.|l|) because pinning-induced order in the liq- 
uid increases with c. The results for B = 3kG are very 
similar, with T c reduced by approximately 3.4K for all 
these values of c. 
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FIG. 1. The phase diagram (T c vs. pin concentration c) 
for B — 2kG. The dashed line is a polynomial fit. Inset: En- 
tropy change As (in units of fcs) and order parameter jump 
Am vs c. Dotted lines are fits to the form const, x (c c — c) 1//2 . 
The approximate position of the critical point is indicated by 
the arrow. 

In Fig^, we have shown the variation of the local den- 
sity p along a line joining two neighboring pinning centers 
for the coexisting crystalline and liquid minima at the 
transition temperature for c = 1/64 and B = 2kG. This 
plot brings out the differences between the structures of 
the coexisting crystal and liquid phases and illustrates 
the ability of our numerical method to provide detailed 
information about the density distribution in inhomoge- 
neous states. In the liquid state, we find [|l9| the expected 
six-fold angular modulation |pj of the density. 

At larger values of c, (/ <6) the behavior found (for 
both values of B) is significantly different. Here, the ap- 
parent minima obtained in "heating" runs (increasing T 
in steps from a crystalline state at low T and using the 
minimum obtained for the last T as the input for the 
current minimization) and in "cooling" runs (decreasing 



T in steps from a liquid state at high T) have almost the 
same free energy, but significantly different density dis- 
tributions [^9| and values of m. At c = 1/25, (I = 5) the 
difference in the values of m peaks at T = T x ~ 21.2K 
for B = 2kG, as shown in Fig.|[ These results suggest 




FIG. 2. The normalized local density p/po along a line 
joining two adjacent pinning centers. Data for the coexist- 
ing crystal (solid line) and liquid (dotted line) for c = 1/64, 
B = 2kG are shown. The peak near x/ao = 1 is at the loca- 
tion of a pinning center. 

that F has only one minimum that becomes very "flat" 
near T = T x : our minimization routine, which assumes 
that a minimum has been reached when the gradient of 
the free energy becomes smaller than a small convergence 
parameter, stops at slightly different points when ap- 
proaching a very flat minimum from different directions. 
To check this, we have calculated the free energy of con- 
figurations {pi(x)} defined by Pi{x) — xp^ + (1 — x)p[ 2 \ 

where {p^ } and {p^} are the configurations obtained 
in heating (crystal) and cooling (liquid) runs, respec- 
tively, and < x < 1 is a mixing parameter. A plot 
of the free energy as a function of x (or of m(x) — 
xm^ 1 ' + (l—x)mP\ where mP'' and m' 2 ' are the order pa- 
rameters in the two configurations) exhibits a minimum 
at x — xo ~ 0.5 at all temperatures. A typical plot, for B 
= 2kG, T = 21. 2K, is shown in the inset of Fig.|. In con- 
trast, similar plots for smaller values of c exhibit a clear 
maximum at an intermediate value of x. These results 
confirm that for c = 1/25 or more, the free energy has a 
unique minimum that lies between the configurations ob- 
tained in heating and cooling runs. Thus, no first-order 
transition occurs at c = 1/25, and the line of first-order 
transitions found for smaller values of c must end at a 
critical point lying between c = 1/36 and c = 1/25. At 
c = 1/25, the change from liquid- like to solid-like behav- 
ior occurs as a sharp crossover. We identify the temper- 
ature T x , which coincides with the temperature at which 
the temperature-derivative of the "equilibrium" value, 



3 



m(xo), of the order parameter peaks, as the crossover 
temperature. The sharpness of the observed crossover 
(see Fig.g) suggests that c = 1/25, T = T x ~ 21. 2K is 
close to the critical point for i? = 2kG, as shown in Fig||. 
For c = 1/16 the crossover is smoother. Our results for 
B = 3kG are very similar, with T x ~ 17. 7K for c = 1/25. 
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FIG. 3. The "heating", "cooling" and "equilibrium" val- 
ues of the order parameter m (see text) as functions of T for 
c = 1/25, B = 2kG. The solid line is a polynomial fit to the 
equilibrium data. Inset: Plot of the free energy vs. m(x) at 
T = 21. 2K. The dotted line is the best fit to Eq.(§). 

We can describe the basic Physics underlying this 
phase diagram as follows. In the presence of commen- 
surate periodic pinning the liquid and the crystal have 
the same symmetry. Since the degree of order in the liq- 
uid increases with c, the liquid and the crystal become 
indistinguishable beyond a critical value of c. Thus, it 
is possible to go from one phase to the other without 
crossing a phase boundary. A simple Landau theory il- 
lustrates this. Symmetry considerations p(J] suggest the 
following Landau expansion for F: 



(3F = —a2m 2 — — a^m 3 + — a^m A 



hm, 
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where 03 and 04 are positive constants, ai is a decreas- 
ing function of temperature, and the "ordering field" h is 
proportional to the pin concentration c. This free energy 
leads to a first-order transition for h < h c — a\/21a\ and 
a critical point at h = h c , = &2c = a§/3a4. The tran- 
sition temperature increases with h, and the latent heat 
and Am vanish as [h c — h) 1 / 2 as h approaches h c from 
below. As shown in Fig.[TJ (inset), our data for As and 
Am are well-described by the form oc (c c — c) 1 / 2 with c c 
close to 1/25. For a more quantitative comparison, we 
have fitted the f3F vs. m data for B = 2kG, T = 21. 2K 
to the form of Eq. (|J) (inset of Fig.||) . The values of a 2c 
and h c , calculated from the best-fit values of 03 and 04, 



are |l9| less than 1.5% lower than the best-fit values of 
ai and h, indicating that the critical point for B = 2kG 
is very close to c = 1/25, T = 21.2K. This explains the 
sharpness of the crossover at c = 1/25. Our results for B 
= 3kG are very similar. The simple Landau theory, thus, 
provides a semi-quantitative account of our density func- 
tional results and suggests that a critical point lies very 
close to c = 1/25 for both values of B. This critical point 
is analogous to the liquid-gas critical point in mean-field 
theory. Fluctuations are expected to change this corre- 
spondence because the symmetry of our order parameter 
is different from that for the liquid-gas transition. 

Theoretical and experimental investigations of the na- 
ture of this critical point in highly anisotropic HTSCs 
would be welcome. The critical point should be experi- 
mentally accessible: the pin lattice spacing for B — 2kG, 
c = 1/25 should be ~ 0.55 /im, close to the spacing of the 
radiation- induced pin array of Ref. |J . Applications of 
our numerical method to other condensed matter systems 
|yl[r2| of a similar nature would also be very interesting. 
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